function p_ri = f_ri(iter,F0,y,X,stations,prob_treat)

    station1 = stations(:,1);
    station2 = stations(:,2);
    station3 = stations(:,3);
    
    sim_data = csvread('stations_randsim.csv',1,0);

    stationsim = sim_data(:,1);
    treatsim = sim_data(:,2:10001);

[N,k] = size(X);

treat1 = zeros(N, iter);
treat2 = zeros(N, iter);
treat3 = zeros(N, iter);

treatany = zeros(N,iter);
wgt=zeros(N,iter);
F = zeros(1,iter);
test = zeros(1,iter);

%Getting ac-wise treatment status
for j = 1:iter
    for i = 1:N
        for s = 1:60 
            if stationsim(s)== station1(i)
                treat1(i,j)=treatsim(s,j);
            end;
            if stationsim(s)==station2(i)
                treat2(i,j)=treatsim(s,j);
            end;
            if stationsim(s)==station3(i)
                treat3(i,j)=treatsim(s,j);         
            end;
            if treat1(i,j)+treat2(i,j)+treat3(i,j)>0
                treatany(i,j)=1;
            end;
        end;
    end;
end;

%Weights
for j = 1:iter
    for i = 1:N
        wgt(i,j)=(treatany(i,j))/prob_treat(i,1) + (1-treatany(i,j))/(1-prob_treat(i,1));
    end;
end;

for a = 1:iter
    [~,~,~,~,F(1,a)] = ipwreg(y,treatany(:,a),X,wgt(:,a),stations,1);
end; 

%p-value
for l = 1:iter
    if F(1,l) > F0
        test(1,l)=1;
    end;
end;
p_ri = mean(test,2);
